function mom = mom_function(data,data_network,fp)
est_data = [];
est_data_network = [];


%Collect simulated data across different neighborhoods and various
%simulations;
for h = 1 : fp.n_rip_sim
    for j = 1:4
        sample = data{j}(:,:,h);
        
        %Because in the last period there is no choice of PS;
        children_under_11 = ( sample(:,fp.ind_data.child_age)<11 );
        est_data = [ est_data ; sample(children_under_11==1,:) ]  ;        
        est_data_network = [ est_data_network ; data_network{j}(:,:,h) ];
        
    end
end

sample_temp =[];

sd_temp_log_skills = nan(fp.n_rip_sim,1);
sd_temp_log_skills_tp1 = nan(fp.n_rip_sim,1);
sd_temp_log_peers = nan(fp.n_rip_sim,1);

for h = 1 : fp.n_rip_sim
    
data_temp = [];

    for j = 1:4
        sample = data{j}(:,:,h);
        %Because in the last period there is no choice of PS;
        sample_temp = [ sample_temp ; sample ]  ;
        data_temp   = [ data_temp   ; sample ]  ;
    end

    sd_temp_log_skills(h,1) = nanstd(log(data_temp(: ,fp.ind_data.child_C) ) );
    sd_temp_log_skills_tp1(h,1) = nanstd( log(data_temp(: ,fp.ind_data.child_C_tp1) ) );
    sd_temp_log_peers(h,1) = nanstd( data_temp(: , fp.ind_data.peers_log)    );
    sd_temp_log_peers_tp1(h,1) = nanstd( data_temp(: , fp.ind_data.peers_log_tp1)    );

end

sd_log_skills = mean(sd_temp_log_skills);
sd_log_skills_tp1 = mean(sd_temp_log_skills_tp1);
sd_log_peers = mean(sd_temp_log_peers);
sd_log_peers_tp1 = mean(sd_temp_log_peers_tp1);

%Construct school-grade fixed effects;
temp_school_grade_fe = [est_data(:,fp.ind_data.n_sim) + max(est_data(:,fp.ind_data.n_sim)).*(est_data(:,fp.ind_data.child_age)-8) ];
temp_school_grade_fe = temp_school_grade_fe - min(temp_school_grade_fe) + 1 ;
temp_fe = dummyvar(temp_school_grade_fe);
school_grade_fe = temp_fe(:,1:end-1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%A) MOMENTS ABOUT INDIRECT INFERENCE ON PS, CHILDREN'S
%SKILLS OR PEERS DYNAMICS.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Indirect Inference on Parenting Style;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y =  est_data( :  ,fp.ind_data.ParStyle);
X=[ ones(size(Y,1),1) log(est_data( : ,fp.ind_data.child_C))./sd_log_skills  (est_data( : ,fp.ind_data.peers_log) - nanmean(est_data( : ,fp.ind_data.peers_log) ) )./sd_log_peers school_grade_fe];

check_var = [Y X ];
check_var = check_var(~any(isinf(check_var),2),:) ;
check_var = check_var(~any(isnan(check_var),2),:) ;

y_var = check_var(:,1);
x_Var = check_var(:,2:end);

temp_coef_PS = ols(y_var ,x_Var );


coef_PS(fp.ind_mom.coef_PS.child_C,1) = temp_coef_PS(2);
coef_PS(fp.ind_mom.coef_PS.peers,1)   = temp_coef_PS(3);
%coef_PS(fp.ind_mom.coef_PS.const,1)   = temp_coef_PS(1);

display('Mean of PS')
display('Model    Data')

mean_PS = nanmean(y_var);

[mean_PS fp.mom_data.mean_PS]

display('OLS Coefficients:  1.child_skills  2.peers_skills 3.constant')
display('Model    Data')

[ coef_PS fp.mom_data.coef_PS]



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Indirect Inference on Dynamics of a Child's Skills;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Y =  ( log(est_data( :  ,fp.ind_data.child_C_tp1)) - nanmean(log(est_data( :  ,fp.ind_data.child_C_tp1))) )./sd_log_skills_tp1 ;


X=[ ones(size(Y,1),1) log(est_data( : ,fp.ind_data.child_C))./sd_log_skills  (est_data( : ,fp.ind_data.peers_log) - nanmean(est_data( : ,fp.ind_data.peers_log) ) )./sd_log_peers est_data( :  ,fp.ind_data.ParStyle) school_grade_fe];

check_var = [Y X ];
check_var = check_var(~any(isinf(check_var),2),:) ;
check_var = check_var(~any(isnan(check_var),2),:) ;

y_var = check_var(:,1);
x_Var = check_var(:,2:end);

temp_coef_skills_tp1 = ols(y_var ,x_Var );

coef_skills_tp1( fp.ind_mom.coef_skills_tp1.child_C ,1) = temp_coef_skills_tp1(2);
coef_skills_tp1( fp.ind_mom.coef_skills_tp1.peers   ,1) = temp_coef_skills_tp1(3);
coef_skills_tp1( fp.ind_mom.coef_skills_tp1.PS      ,1) = temp_coef_skills_tp1(4);
%coef_skills_tp1( fp.ind_mom.coef_skills_tp1.const   ,1) = temp_coef_skills_tp1(1);


display('Next Period Skills : 1.child_skills  2.peers_skills 3.Parenting_Style ')
display('Model    Data')
[ coef_skills_tp1 fp.mom_data.coef_skills_tp1]



y_var_ps1 = y_var( x_Var(:,4)==1 );
x_Var_ps1 = x_Var( x_Var(:,4)==1 ,  [ 1 2 3 5:end ]);
temp_coef_skills_tp1_ps1 = ols(y_var_ps1 ,x_Var_ps1 ) ;


y_var_ps0 = y_var( x_Var(:,4)==0 );
x_Var_ps0 = x_Var( x_Var(:,4)==0 ,  [ 1 2 3 5:end ]);
temp_coef_skills_tp1_ps0 = ols(y_var_ps0 ,x_Var_ps0 ) ;

coef_skills_tp1_ps0( fp.ind_mom.coef_peers_tp1.child_C ,1)= temp_coef_skills_tp1_ps0(2);
coef_skills_tp1_ps0( fp.ind_mom.coef_peers_tp1.peers   ,1)= temp_coef_skills_tp1_ps0(3);

coef_skills_tp1_ps1( fp.ind_mom.coef_peers_tp1.child_C ,1)= temp_coef_skills_tp1_ps1(2);
coef_skills_tp1_ps1( fp.ind_mom.coef_peers_tp1.peers   ,1)= temp_coef_skills_tp1_ps1(3);


display('Next Period Skills (PS=0) : 1.child_skills  2.peers_skills    ')
display('Model    Data')
 [ coef_skills_tp1_ps0 fp.mom_data.coef_skills_tp1_ps0  ]
 
display('Next Period Skills (PS=1) : 1.child_skills  2.peers_skills    ')
display('Model    Data')
 [ coef_skills_tp1_ps1 fp.mom_data.coef_skills_tp1_ps1  ]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Indirect Inference on Dynamics of Peer's Quality;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y = ( est_data( :  ,fp.ind_data.peers_log_tp1) - nanmean(est_data( :  ,fp.ind_data.peers_log_tp1)) )./sd_log_peers_tp1 ;
X=[ ones(size(Y,1),1) log(est_data( : ,fp.ind_data.child_C))./sd_log_skills  (est_data( : ,fp.ind_data.peers_log) - nanmean(est_data( : ,fp.ind_data.peers_log) ) )./sd_log_peers est_data( :  ,fp.ind_data.ParStyle)  school_grade_fe];

check_var = [Y X ];
check_var = check_var(~any(isinf(check_var),2),:) ;
check_var = check_var(~any(isnan(check_var),2),:) ;

y_var = check_var(:,1);
x_Var = check_var(:,2:end);

temp_coef_peers_tp1 = ols(y_var ,x_Var );

coef_peers_tp1( fp.ind_mom.coef_peers_tp1.child_C ,1)= temp_coef_peers_tp1(2);
coef_peers_tp1( fp.ind_mom.coef_peers_tp1.peers   ,1)= temp_coef_peers_tp1(3);
coef_peers_tp1( fp.ind_mom.coef_peers_tp1.PS      ,1)= temp_coef_peers_tp1(4);
%coef_peers_tp1( fp.ind_mom.coef_peers_tp1.const   ,1)= temp_coef_peers_tp1(1);


display('Next Period Peers : 1.child_skills  2.peers_skills   3.Parenting_Style 4.constant  ')
display('Model    Data')
[ coef_peers_tp1 fp.mom_data.coef_peers_tp1 ]

y_var_ps1 = y_var( x_Var(:,4)==1 );
x_Var_ps1 = x_Var( x_Var(:,4)==1 ,  [ 1 2 3 5:end ]);
temp_coef_peers_tp1_ps1 = ols(y_var_ps1 ,x_Var_ps1 ) ;


y_var_ps0 = y_var( x_Var(:,4)==0 );
x_Var_ps0 = x_Var( x_Var(:,4)==0 ,  [ 1 2 3 5:end ]);
temp_coef_peers_tp1_ps0 = ols(y_var_ps0 ,x_Var_ps0 ) ;

coef_peers_tp1_ps1( fp.ind_mom.coef_peers_tp1.child_C ,1)= temp_coef_peers_tp1_ps1(2);
coef_peers_tp1_ps1( fp.ind_mom.coef_peers_tp1.peers   ,1)= temp_coef_peers_tp1_ps1(3);

coef_peers_tp1_ps0( fp.ind_mom.coef_peers_tp1.child_C ,1)= temp_coef_peers_tp1_ps0(2);
coef_peers_tp1_ps0( fp.ind_mom.coef_peers_tp1.peers   ,1)= temp_coef_peers_tp1_ps0(3);

display('Next Period Peers (PS=0) : 1.child_skills  2.peers_skills    ')
display('Model    Data')
[ coef_peers_tp1_ps0 fp.mom_data.coef_peers_tp1_ps0 ]

display('Next Period Peers (PS=1) : 1.child_skills  2.peers_skills    ')
display('Model    Data')
[ coef_peers_tp1_ps1 fp.mom_data.coef_peers_tp1_ps1 ]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Indirect Inference on Investments;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sample_PS0 = ( est_data( :  ,fp.ind_data.ParStyle)==0 );
sample_PS1 = sample_PS0==0; 

Y = est_data(:, fp.ind_data.inv) ;
Y = (Y - nanmean(Y) )  ;

X=[ ones(size(Y,1),1) log(est_data( : ,fp.ind_data.child_C))./sd_log_skills  (est_data( : ,fp.ind_data.peers_log) - nanmean(est_data( : ,fp.ind_data.peers_log) ) )./sd_log_peers   school_grade_fe];

%If PS==0;
check_var = [Y X ];
check_var = check_var( sample_PS0==1 , :);
check_var = check_var(~any(isinf(check_var),2),:) ;
check_var = check_var(~any(isnan(check_var),2),:) ;

y_var = check_var(:,1);
x_Var = check_var(:,2:end);

temp_coef_inv = ols(y_var ,x_Var );

coef_inv_PS0( fp.ind_mom.coef_inv.child_C , 1) = temp_coef_inv(2);
coef_inv_PS0( fp.ind_mom.coef_inv.peers   , 1) = temp_coef_inv(3);
%coef_inv_PS0( fp.ind_mom.coef_inv.const   , 1) = temp_coef_inv(1);


display('Parental Investments (PS=0) : 1.child_skills  2.peers_skills   3.Parenting_Style 4.constant  ')
display('Model    Data')
[ coef_inv_PS0 fp.mom_data.coef_inv_PS0]

%If PS==1;
check_var = [Y X ];
check_var = check_var( sample_PS1==1 , :);
check_var = check_var(~any(isinf(check_var),2),:) ;
check_var = check_var(~any(isnan(check_var),2),:) ;

y_var = check_var(:,1);
x_Var = check_var(:,2:end);

temp_coef_inv = ols(y_var ,x_Var );

coef_inv_PS1( fp.ind_mom.coef_inv.child_C , 1) = temp_coef_inv(2);
coef_inv_PS1( fp.ind_mom.coef_inv.peers   , 1) = temp_coef_inv(3);
%coef_inv_PS1( fp.ind_mom.coef_inv.const   , 1) = temp_coef_inv(1);


display('Parental Investments (PS=1) : 1.child_skills  2.peers_skills   3.Parenting_Style 4.constant  ')
display('Model    Data')
[ coef_inv_PS1 fp.mom_data.coef_inv_PS1]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Mean of Child Skills by Age;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

mean_skills(1,1) = nanmean( log(sample_temp( sample_temp( : ,fp.ind_data.child_age,1)==9  ,fp.ind_data.child_C,1)) ) ;
mean_skills(2,1) = nanmean( log(sample_temp( sample_temp( : ,fp.ind_data.child_age,1)==10 ,fp.ind_data.child_C,1)) ) ;
mean_skills(3,1) = nanmean( log(sample_temp( sample_temp( : ,fp.ind_data.child_age,1)==11 ,fp.ind_data.child_C,1)) ) ;
mean_skills(4,1) = nanmean( log(sample_temp( sample_temp( : ,fp.ind_data.child_age,1)==11 ,fp.ind_data.child_C_tp1,1)) ) ;
display('Mean Child Skills by Age')
display('Model    Data')
[mean_skills  fp.mom_data.mean_skills ]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Mean of Investments (by PS);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

sample_PS0 = ( est_data( :  ,fp.ind_data.ParStyle)==0 );
sample_PS1 = sample_PS0==0; 

Y = est_data(:, fp.ind_data.inv) ;
Y = (Y - nanmean(Y) )  ;

mean_inv(1,1)=nanmean(Y(sample_PS0==1));
mean_inv(2,1)=nanmean(Y(sample_PS1==1));

display('Mean Investments by PS')
display('Model    Data')
[mean_inv  fp.mom_data.mean_inv ]

%%%%%%%%%%%%%%%%%;
% N of Friends;
%%%%%%%%%%%%%%%%%;

data_n_friends = est_data_network( : , fp.ind_data_network.n_friends );
mean_nfriends = nanmean(data_n_friends(data_n_friends>0)) ; 

display('Average N of Friends')
display('Model    Data')
[mean_nfriends  fp.mom_data.mean_nfriends ]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Longitudinal Regressions on PS;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

delta_skills = [];
delta_peers  = [];
delta_PS     = [];
delta_inv    = [];
author_tm1   = [];

for g = 10:1:11
    g_m1 = g-1;
    
        children_age_t = ( sample_temp(:,fp.ind_data.child_age)==g );
        children_age_t_m1 = ( sample_temp(:,fp.ind_data.child_age)==g_m1 );    
                       
        delta_skills = [ delta_skills ;  ( log(sample_temp( children_age_t==1 ,fp.ind_data.child_C)) - log(sample_temp( children_age_t_m1==1 ,fp.ind_data.child_C)))./ nanstd(log(sample_temp( : ,fp.ind_data.child_C)))  ];
        delta_peers  = [ delta_peers ;  (  sample_temp( children_age_t==1 ,fp.ind_data.peers_log)  -  sample_temp( children_age_t_m1==1 ,fp.ind_data.peers_log) )./nanstd(sample_temp( : ,fp.ind_data.peers_log))  ];
        delta_PS     = [ delta_PS   ; (sample_temp( children_age_t==1   ,fp.ind_data.ParStyle) - sample_temp( children_age_t_m1==1 ,fp.ind_data.ParStyle) ) ]  ; 
        delta_inv    = [ delta_inv  ; (sample_temp( children_age_t==1, fp.ind_data.inv) - sample_temp( children_age_t_m1==1, fp.ind_data.inv) ) ] ; 
        author_tm1   = [ author_tm1 ;  sample_temp( children_age_t_m1==1 , fp.ind_data.ParStyle ) ] ;                       

end


Y = delta_PS;
X = [ ones(size(delta_skills)) delta_skills delta_peers ];
check_var = [Y X ];
check_var = check_var(~any(isinf(check_var),2),:) ;
check_var = check_var(~any(isnan(check_var),2),:) ;

y_var = check_var(:,1);
x_Var = check_var(:,2:end);

temp_coef_PS_longit = ols(y_var ,x_Var );

coef_PS_longit(fp.ind_mom.coef_PS_longit.child_C,1) = temp_coef_PS_longit(2);
coef_PS_longit(fp.ind_mom.coef_PS_longit.peers,1) = temp_coef_PS_longit(3);

display('Delta PS: : 1.Delta child_skills  2.Delta peers_skills')
display('Model    Data')
[coef_PS_longit  fp.mom_data.PS_longit ]



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Longitudinal Regressions on Investments;
%(Conditional on being authoritarian at (t-1) )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

Y = delta_inv;
X = [ ones(size(author_tm1))  author_tm1  ];
check_var = [Y X ];
check_var = check_var(~any(isinf(check_var),2),:) ;
check_var = check_var(~any(isnan(check_var),2),:) ;

y_var = check_var(:,1);
x_Var = check_var(:,2:end);

temp_coef_inv_longit = ols(y_var ,x_Var );

display('Delta Inv: : 1.Constant  2.Authoritarian (t-1)')
display('Model    Data')
[temp_coef_inv_longit(1:2)  fp.mom_data.Inv_longit  ]


mean_PS_neighborhoods = NaN(4,1);
ind = 1;
for j = 1:1:4
   
    mean_PS_neighborhoods(ind,1) = nanmean(sample_temp(sample_temp(:,fp.ind_data.neighborhood)==j , fp.ind_data.ParStyle )); 
    
    ind = ind + 1 ;
    
end


display('Mean PS by Neighborhood')
display('Model    Data')
[mean_PS_neighborhoods  fp.mom_data.mean_PS_between  ]


sim_moments=  [ mean_PS;  coef_PS; coef_skills_tp1 ;coef_skills_tp1_ps0 ; coef_skills_tp1_ps1 ; coef_peers_tp1; coef_peers_tp1_ps0 ; coef_peers_tp1_ps1 ; coef_inv_PS0; coef_inv_PS1 ; mean_skills; mean_inv; mean_nfriends ; coef_PS_longit ; temp_coef_inv_longit(1:2) ; mean_PS_neighborhoods ] ; 

mom = sim_moments;
    